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Abstract. - A new method is used to measure the stress and elastic constants of purely 

entropic phantom networks, in which a fraction p of neighbors are tethered by inextensible 

. ■ bonds. We find that close to the percolation threshold pc the shear modulus behaves as (p—pc), 

where the exponent / « 1.35 in two dimensions, and / « 1.95 in three dimensions, close 

e, to the corresponding values of the conductivity exponent in random resistor networks. The 

components of the stiffness tensor (elastic constants) of the spanning cluster follow a power law 
i-p^ ' ^ (p — pcY , with an exponent g ~ 2.0 and 2.6 in two and three dimensions, respectively. 

o 

In the gelation process, monomers or short polymers in a fluid solution are randomly cross- 
linked. At a certain moment during the reaction, a macroscopically large network, the gel, 
spans the system. At this point, the system changes from a fluid-like (sol) to a solid- like (gel) 
phase that has a finite shear modulus. The geometry of gels is frequently described by the 
percolation model 13] . The percolation geometry is usually defined on a lattice, by randomly 
^\ . occupying a fraction p of the bonds (or sites). The gel point is identified with the percolation 
^^ ' threshold pc, the critical bond (site) concentration above which a spanning cluster is formed. 
^^ . Percolation theory predicts that close to Pc quantities like the average cluster size or the gel 
fraction have power laws dependence on (p — pc) with universal exponents, some of which have 
t^ . been measured experimentally for gel systems [g|. 

Near the sol-gel transition typical polymer clusters are very large, tenuous and floppy. 
Elastic properties of such systems are primarily determined by the entropy, i.e., distortions 
'^ ' of a sample barely modify its energy, but they decrease the available phase space (decrease 
C , entropy) and, thus, increase the free energy. Like geometrical quantities near pc, the shear 
modulus is also expected to follow a power law: fi ~ (p — pc)^. De Gennes [|| used an analogy 
between gel elasticity and conductivity of random resistor networks (RRN), and conjectured 
that the exponent / should be equal to the exponent t describing the conductivity S of RRN 
S^ . close to pc^ ^ ^ {p ~ PcY ■ Alternative theories take different approaches and lead to different 
$H ' exponents [||. An exact calculation of the critical behavior of p, which takes into account 
excluded volume (EV) and entanglements effects, is not yet available. Experimental values of 
/, measured for different polymeric systems, are very scattered [^. One of the reasons for the 
variety of the experimental results is the mixing of the entropic and energetic contributions 
to the gel elasticity, which influences the "effective" exponent. 
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Neglect of EV interactions, i.e., treating a phantom system, may strongly modify the 
physics. Nevertheless, it is frequently done either because in certain situations (such as dense 
polymer melts) EV interactions effectively cancel out B, or because from purely theoreti- 
cal point of view phantom systems are more tractable and may serve as a starting point for 
studying real systems. Phantom systems maintain the correct connectivity, which is one of 
the important characteristics of a polymer network. A feature common to most phantom 
networks (independently of the detailed shape of the microscopic potential) is the fact at zero 
tension the probability density that two distant nodes are separated by r, takes a Gaussian 
form ~ exp [—5^''^] • For linear polymers this is a consequence of the central limit theorem, 
while for more complicated systems this can be demonstrated numerically m. The thermo- 
dynamic behavior of a phantom network can, therefore, be described very accurately within 
a phantom Gaussian network (PGN) model, in which each bond of the network is replaced 
by a Gaussian spring having the energy E = ^Kr^, where r is its end-to-end distance |8|. 
Corrections to Gaussian behavior can be studied by considering a phantom nearly-Gaussian 
network (PNGN), in which the springs' energies include an additional small term equal to 
jar'^. In this paper we describe results of a numerical study of the elasticity of tethered 
phantom percolating networks, and compare our results with the two models [g[. We use a 
recently developed formalism ||l^ which enables direct calculation of entropic elastic constants 
of tethered systems. We show that the shear modulus behaves near pc like the conductivity of 
RRN as predicted by the PGN model, while the elastic stiffness tensor of the spanning cluster, 
which according to PGN model is supposed to vanish Q, also exhibits a power law behavior 
near pc with a significantly larger critical exponent. The last result is a consequence of the 
deviation from the Gaussian behavior, and can be understood within the PNGN model Q. 

In a homogeneous deformation we distort the boundaries of a system is such a way that a 
separation R between a pair of surface points is modified into r, which is linearly related to 
R via a position-independent matrix. In this case the new squared distance r^ = RiRj{Sij + 
'^rjij), where the subscripts denote Cartesian coordinates, dij is the Kronecker delta, rjij is 
the Lagrangian strain tensor, and summation over repeated indices is implied. The stress 
and elastic stiffness tensors, ct^ and Cijki , respectively, are then defined as the coefficients of 
the expansion of the free energy density in the strain variables: /({jy}) = /({O}) + o'ijrjij + 
^CijkirjijTjki + ■ ■ ■ Close to Pc, percolation networks "forget" the details of the lattice and 
behave like isotropic systems, and therefore the stress tensor a^x = '^yy — <^zz = cr = — P, 
where P is the pressure. Isotropic systems have only three different non-vanishing elastic 

constants. On ^ ^xxxx — ^yyyy — ^zzzz i ^12 ^ ^xxyy — ^yyzz — ^zzxx — ■ ■ • 5 aUQ 

C44 = ^{Cxyxy + Cxyyx) = \{Cyzyz + Cyzzy) = ■■■■, which are related by: Cn = C12 + 2C44. 
Frequently, one finds it more useful to describe the elastic behavior in terms of the shear 
modulus /x = C44 — P, and the hulk modulus n = ^(Cn -f C12) [for two-dimensional (2D) 
systems], or k = |(Cii -(- 2Ci2 + P) [for three-dimensional (3D) systems]. In a percolation 
phantom system the contribution of the different clusters are additive. Each finite cluster, not 
connected to the boundaries of the system, contributes as a single atom of an ideal gas. Thus 
A^o free finite clusters confined within volume V at temperature T, produce stress equal to 

a = ^^, and elastic constant C44 = ^'^'^ , where k is the Boltzmann constant. Although 

both C44 and P are affected by the presence of finite clusters, we observe that the (ideal 
gas) contribution of the finite clusters cancels out in the definition definition of /i. Since 
finite clusters play such unremarkable role in the problem of elasticity, we will disregard them 
completely, and in the remainder of this work the stress, elastic constants and elastic moduli 
will refer to the contribution of the spanning cluster alone. The latter depends on the details 
of the potential and connectivity. However, in the case of PGN few simple properties exist [^: 
(1) Since the energy of a Gaussian spring is proportional to r^, while the squared distances 
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Fig. 1 - Part of initial (a) and equilibrated (b) configurations of the 2D system (p — 0.405, b = 1.05). 



between the points on the boundaries are linear in 77^, it can be shown the free energy of 
the spanning cluster does not include quadratic terms in 77^, and therefore its Cijki = 0. 
(2) The stress tensor is equal to the conductivity tensor of an equivalent resistor network in 
which each spring of the spanning cluster with a force constant K is replaced by a resistor of 
conductance K. For isotropic PGNs jj, — a — 'S, where S is the conductivity of the equivalent 
resistor network, and we, thus, find that f — t. In the PNGN model near pc, a and /j, are 
still dominated by the Gaussian term rather than by the non-Gaussian perturbation, and 
we recover the equality f — t. Non-Gaussian corrections are manifested by non-vanishing 
elastic constants, which for percolation PNGNs are expected to behave as C ~ (p — Pc)^ with 

9>fM- 

In this work we investigate the elastic behavior of networks for which the problem of 
mixing of the entropic and energetic components does not exist. Our systems consists of 
point-like atoms connected by "tethers" that have no energy, but simply limit the distance of 
a connected pair to be smaller than some value b. Since the internal (potential) energy of the 
system vanishes, its thermodynamic behavior is purely entropic. We generated the (quenched) 
topologies by considering bond percolation problem on 2D triangular {pc = S ~ 0.349) and 
3D faced-centered-cubic (pc — 0.12) lattices, with a fraction p of bond present. Each present 
bond was replaced by a tether, while each site became an "atom" without EV, and the system 
was allowed to move in continuum. Fig. m (a) depicts an initial 2D configuration of the system, 
which equilibrates into configuration of the kind depicted in fig. |] (b). As expected, finite 
clusters and dangling ends of the spanning cluster contract relative to their linear size in the 
initial quenched construction [O. The size of the backbone, on the other hand, is fixed by 
the boundary conditions and, therefore, it looks like a collection of loops of the size of the 
percolation correlation length. 

The non-interacting character of phantom networks significantly simplifies the numerical 
procedure: (1) Since we are not interested in the (trivial) contribution of the finite clusters 
they were removed from the simulations. (2) Dangling ends of the spanning cluster do not 
contribute neither to the stress nor to the elastic constants and, therefore, they can also be 
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Fig. 2 - Logarithmic plot of the stress a and the shear modulus ^ as a function of (p — pc), for 2D 
systems. The slope of the solid line is 1.35. Results are in kT/b^ units. 

Fig. 3 - Logarithmic plot of the elastic constants Cn and C44 as a function of (p — pc), for 2D systems. 
The slope of the lines is « 2. Results are in kT/b^ units. 



removed. Thus, for every quench we identified the backbone (using the "burning" algorithm 
|I3|, which was slightly modified to deal with the periodic boundary conditions applied in 
the simulations), and explored the configuration space using a Monte Carlo (MC) updating 
scheme |0| in which the conventional Metropolis single atom steps are replaced by collective 
steps of chains of atoms. At each MC time unit we made a number of move attempts (with 
acceptance probability ~ 0.5) equal to the number of atoms. In the 2D simulations, we used 
a 120 X 138 triangular lattice (that has an aspect ratio very close to 1) with nearest-neighbor 
spacing bo = 1, and a number of quenched topologies that ranged from Nt = 200 for p closest 
to Pc, down to Nt = 20 far from pc- In the 3D simulations we used systems of 24^ cubic 
unit cells (each containing 4 atoms), i.e., of linear size L = 24-\/2&o, with nearest-neighbor 
spacing bo = 1, and 30 < Nt < 150. The duration of the MC run of each individual sample 
was at least 50 times larger than the relaxation time which we estimated from the expression 
r = dkTL^p/{TT^^s^), where s is the (average) distance an atom moves in one MC time unit, p 
is the number density of atoms, and d is the dimensionality of the system [[l4| . The value of fi in 
this expression was taken, a posteriori, from the simulations. We used a new method enabling 
the direct measurement of the stress and elastic constants from the probability densities of 
finding maximally extended tethers 1 10 . The error estimates are affected by the fluctuations 



in the values of the measured quantities between the different quenches, and to lesser extent 
by the thermal uncertainties within each sample. The error bars appearing in the graphs 
correspond to one standard deviation of the average. 

In our system we can vary only two non-trivial parameters: bond concentration p, and the 
maximal tether length b (measured in the units of the nearest-neighbor spacing bo). Fig. || 
depicts our results for a and /z as a function of {p — pc) for 2D systems with b = 1.05 and 
b = 1.25. It clearly demonstrates that close to Pc, the network becomes Gaussian: First, 
the difference between p and a decreases as we approach pc, which implies that the elastic 
constant C44 = p — a vanishes faster than both quantities. Second, when plotted in kT/b^ 
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Fig. 4 - Logarithmic plot of the stress a and the shear modulus /x as a function of ( 
systems with b = 1.05. The slope of the solid line is 1.95. Results are in kT/b^ units. 



- Pc), for 3D 



Fig. 5 - Logarithmic plot of the elastic constants Cn, C12, and C44 as a function of {p — pc), for 3D 
systems with b = 1.05. The slope of different lines is ~ 2.65. Results are in kT/b"^ units. 



units, the values of a and /i in systems with different b converge towards each other. This is 
explained by the facts that (a) the stress of a 2D PGN depends only on the topology of the 
network and the value of the springs force constant K; and (b) for the tethered networks, the 
effective K is proportional to kT/})^ [|j. Third, the value of / extracted from the the graphs 
is / = 1.35±0.10, very close to the value of the conductivity exponent t — 1.297±0.007 in 2D 
|15[ . Similar result for the exponent / has been obtained by Plischke et al. |jlra. They used 
central force networks in which both entropy and energy contribute to the elastic properties 
and, by examining systems at several temperatures, removed the energetic component. Close 
to Pc elasticity of central force systems is completely dominated by entropy, and their result 
for / reflects this fact. 

In fig. ^we present our results for the elastic constants Cn and C44, which are supposed to 
vanish in the purely Gaussian case. Shorter tethers correspond to larger values of the elastic 
constants, since they represent more stretched networks, which exhibit stronger deviations 
from Gaussian behavior. Despite almost an order-of-magnitude difference between Cn and 
C44 for the same h, and half an order-of-magnitude difference between the same constants for 
the different values of 6, all the results can be described by a power law [p — PcY , with the 
same exponent g = 2.0 ± 0.2, which is significantly larger than /. (We do not show the elastic 
constant C12, which has large statistical uncertainties that prevent exact determination of the 
power law. The results are, however, consistent with the power laws for the other constants.) 
To further ascertain the universality of g, one would need to increase b to even larger values. 
This, however, would further decrease the values of the elastic constants which in our method 
of simulations |10| would increase the statistical uncertainties, and require increase of the 
simulation length beyond our computational ability. 

Our results for the 3D networks with b = 1.05 are shown in figures ^ and ^ Again, the 
validity of the PGN model is supported by the observation that a and /x converge towards 
each other as we approach pc, following power laws with / = 1.95 ± 0.05, which agrees with 
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the conductivity exponent t ~ 2.003 ± 0.047 in 3D [O. The elastic constants also follow 
power laws with an exponent g = 2.65 ± 0.15. Note that our results confirm the relation 
Cii = Ci2 + 2C44, which indicates that close to Pc, percolating networks behave as isotropic 
systems. At p = I the system has a lower (cubic) symmetry, and there is a gradual deviation 
from this relation with increasing p beyond the regime shown in fig. ra. 

In ref. H we used PNGN model to derive bounds on the exponent g, which is a consequence 
of non-Gaussian behavior: 3t — 2v{d— 1) < g < 4(t— 1) — v{'id — A), where v is the correlation 
length exponent. The perturbative derivation of the bounds was self-consistent only for d > 3. 
We note that our result for 5 in d = 3 is, within the statistical uncertainty, consistent with 
the bounds 2.48 < g < 2.60. One should keep in mind, however, that the PNGN model 
assumes that the coefficient of the quartic perturbation term to the Gaussian spring energy 
is a constant number, while for the tethered network model its effective value may depend on 
the mean stress and, thus, on the position of the bond in the network. 

In conclusion, we studied the critical elastic behavior of purely entropic phantom model 
with topology of a percolating network. The microscopic tethering potential is very different 
from a Gaussian spring. Nevertheless, diluted networks become very "floppy" so that the 
potentials become effectively Gaussian and, consequently, the shear modulus behaves as the 
conductivity of RRN. Non-trivial power law dependence of the elastic constants of the spanning 
cluster on (p—pc) is a signature of a deviation from the Gaussian behavior, and is controlled by 
a critical exponent significantly larger than the exponent of conductivity. Since g characterizes 
a "sub-leading" behavior, a detailed study of a broad class of potentials is needed to verify its 
universality. 
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